%% DoubleAngleB1Map Function: Calculates B1 map from double angle sequence data.
% 
% Function inputs two images with a flip angle ratio of alpha and 2*alpha. 
% The nominal flip angle is a required input, and is used to make the B1 
% map values a relative, to be used for image calibration.
%
% Mathematical Proof
% 
% From geometrical arguments:
% Image1 = Mz0*sin(alpha), Image2 = Mz0*sin(2*alpha)
%
% Thus, Mz0 = Image1/sin(alpha) = Image2/sin(2*alpha)
% 
% sin(2*alpha)/sin(alpha) = Image2/Image1 where this ratio can has a 
% maximum of 2 for very small angles (or close to n*2*pi).
%
% Trigonometric Identity -> sin(2*theta)=2*sin(theta)*cos(theta)
%
% Thus,
% 
% Image2/Image1 = 2*sin(alpha)*cos(alpha)/sin(alpha)
%
% and
%
% alpha = acos(1/2*Image2/Image1) Q.E.D.

%% Function Declaration
%
function [ b1map ] = DoubleAngleB1Map( doubang1, doubang2, NominalAngle)
%% Cast int32 images to double.
% Needed due to use of loadminc tool.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
doubang1=cast(doubang1, 'double');
doubang2=cast(doubang2, 'double');
%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate the cosine of the true flip angles.
%
r = abs(doubang2./doubang1);      

cos_arg = 1/2*r;

%% Filter out cases where r > 2.
% r should/can not be larger than than two, except in the case of noise.
cos_arg = double(cos_arg).*(r<=2) + ones(size(r)).*(r>2);

%% Calculate true flip angle.
%
alpha = acos(cos_arg); %alpha is in radians

%%  Convert flip angle map to relative B1 map.
% To convert alpha to a relative factor, we divide by the nominal flip 
% angle in radians.

b1map = alpha/(NominalAngle*pi/180);

end

